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Abstract 


This thesis deals with the topic of generalized finite element method (GFEM) with particular 
reference to linear elastic fracture mechanics (LEFM) The presence of a crack in the domain 
leads to singular stresses at the crack tip, for linear elastic analysis In this case, the con\en- 
tional finite element analysis is known to lead to inferior quality crack parameters The quaht\ 
of extracted crack parameters can be improved only by significant mesh refinements at the crack 
tip Use of refined meshes increases the computational load In the present study, an effort 
has been made to characterise the crack parameters more accurately and (computationally) effi- 
ciently The basic idea is to use specialized basis functions at the crack tip, using the asymptotic 
crack solutions The method uses the partition of unity approach to make the approximation 
conforming The code has been benchmarked with a large set of sample problems Further, a 
crack propagation problem in planar elasticity is considered Lastly, the problem of propagation 
of multiple, arbitrarily placed, cracks is considered 
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Chapter 1 
Introduction 


1.1 Introduction 

The generalized finite element method (GFEM) is a combination of classical finite element 
method (FEM) and partition of unity method (PUM) In generalized finite element (GFEM), 
the standard finite element spaces are augmented by adding special functions which reflect the 
known information about the boundary value problem and input data, e g the singular solution 
obtained from the local asymptotic expansion of the exact solution in the neighbourhood of a 
centre point, etc The essential idea in this method is to add enrichment functions to the finite 
element approximation which contains a discontinuous displacement field The method exploits 
the partition of unity property of finite elements which was noted by Melenk and Babuska [1], 
and Duarte and Oden [2], namely that the sum of the shape functions must be unity The special 
functions are multiplied with the partition of unity to construct an augmented conforming finite 
element space In this way the local approximability afforded by the special functions is included 
in the standard fini te element approximation, while maintaining the existing infrastructure of 
finite element codes 

1.2 Literature survey 

In the field of computational mechanics, extensive research has been done m the effective use of 
the finite element method for problems m fracture mechanics Effort has been primarily focussed 
on the modification of the approximation, incorporating the singularity present in the domain 
Of the various elements developed, collapsed quadrilateral elements and quarter point element 
requires special mention [16] While direct application of standard elements is the most straight- 
forward, a high degree of mesh refinement is required near the crack tip to capture the singular 
stress fields 

Recently, various methods that attempt to do away with the mesh have become popular for 
solving boundary value problems, for example, the element-free galerkin methods of Belytschko 
et al [11], the hp Clouds method of Duarte and Oden [2] These methods in their purest form 
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get rid of the finite element method (FEM) completely, and intend to replace traditional finite 
element codes with a new infrastructure. The methods have had limited success because they 
have not successfully addressed the problem of the numerical integration of the entries of the 
discrete equations. For example, all meshless methods must still employ an integration mesh for 
the numerical integration. However, because the regions of integration which are often intersect- 
ing disks in two-dimensions (2D) or spheres in three-dimensions (3D) do not coincide with the 
integration mesh, it is very difficult to control the error of the numerical integration. Further, the 
various meshless methods also require special techniques for applying essential boundary condi- 
tions, e.g., the method of Lagrange multipliers which raise the issues of stability. An enriched 
element free galerkin method (EFG) was proposed by Fleming et al [4]. However, it has yielded 
good results only when adequate refinement (in the background mesh) is used near the crack tip. 
Moreover, the stresses computed also tend to be oscillatory near the crack tip unless substantial 
refinement is used. All these modern meshless methods are based implicitly or explicitly on using 
a partition of unity over the the domain to ensure continuity of the approximation. The partition 
of unity method (PUM) was developed by Babuska et al [1]. 

Recently, the generalized finite element method (GFEM) was proposed by Strouboulis et al [6] 
to do away with the problems encountered in common meshless methods. The Computational 
Mechanics Company (COMCO) is currently developing a new fast and robust solver for struc- 
tural mechanics problem. It is based on a new solver technology, which promises to combine the 
modeling advantages and flexibility of finite element methods with the ease of use of meshless 
methods and with computational performance compatible or faster than modern finite element 
codes. The solver is based on proprietary extensions of meshless cloud methods and generalized 
finite element method (GFEM). It is being implemented in an hp-adaptive context within the 
computational environment PHLEX. Formally, the method is named generalized finite element 
partition of unity method or, for short, an element partition method (EPM) [7]. 

The present work, is a first step towards developing a generalized finite element based analysis 
platform for the study of quasi-static and dynamic crack propagation. The local asymptotic 
solution, at the crack tip, is chosen for enrichment. These special functions are made conforming 
by employing partition of unity based cut-off functions. 

1.3 Thesis organization 

The thesis is organised as follows: 

1. In chapter 2, we briefly review the concepts of the re-entrant corner and its generalisation 
to linear elastic fracture mechanics. The procedure employed to calculate the different 
crack parameters are outlined in this chapter. 

2. In chapter 3, a general procedure is implemented, based on energy principle, for deriving 
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the weak formulation studied in this work 


3 In chapter 4, we discuss the mam features of the generalized finite element method employed 
in this work An adaptive integration scheme tailored to the present work has been detailed 
elaborately 

4 In chapter 5, we present the numerical results obtained and a qualitative analysis has been 
done to compare the classical finite element solution and generalized finite element solution 
Lastly an effort is made to predict the crack path using the above principle 

5 Finally m chapter 6, we summarize our work and draw the relevant conclusions and the 
scope for future work 


/ 
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Chapter 2 

Linear elastic fracture mechanics 

2.1 Two dimensional planar elasticity 

Equilibrium Equation 

In planar elasticity problems, the variation of stress components in third direction is negligible 
Also the body force part in fracture mechanics is absent, hence the equilibrium equations can be 
simplified as follows 


dan dai2 
dxi dx2 
da 12 da 22 
dxi dx2 


0 

0 


(2 1 ) 
(2 2 ) 


Strain displacement relationship 


If ui and U2 are the two planar displacement components, then the strain-displacement relation 


IS given by 


£■11 

£22 

£12 


dui 


dx\ 

du2 

dx2 

1 /dui 

2 \dx2 


+ 



The compatibility relation is given by 

d'^Sn d'^S22 _ ^ d'^en 
dxi dxj dxidx2 


0 


(2 3 ) 
(2 4 ) 
(2 5 ) 

(2 6 ) 


Stress strain relationship 


For the linear elasticity problem with isotropic materials, the stress strain relationship can be 


given as 


£11 = ~ ^(<^22 + <733)] 


(2 7 ) 
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£22 = 

^[(722 - 1^(C^11 

+ ^^ 33 )] 

(2 8) 

£33 = 

^[<733 - i^(o-ii 

+ 0'22)] 

(2 9) 


(1 + I/) 

E = 

1 

(2 10) 

£12 = 



Plane deformation 

For the case of plane stress i e for a very thin plate, the out of plane stresses are negligible in 
the interior points of the plate Thus we have 


(Ji3 — a23 = 0‘33 = 0 

Thus the stress strain relation reduces to 

Sii = -^[(Tu ~ ^<^22] 

£22 = -j^[(^22 — ^CTii] 

1 

£12 - 2511.2 


(2 11 ) 

(2 12 ) 
(2 13) 
(2 14) 


Similarly for the plane strain case, 1 e for sufficiently thick plates , the through the thickness 
strain leads to 

£^13 = £23 = £33 = 0 (2 15) 

(733 = + <^ 22 ) (2 16) 

Thus the stress strain relation reduces to 


- (^<^221 


(1 - 


£22 = l^[.22-j^...l 

or, a general representation for the planar constitutive relationship can be given by 

1 r * 1 


■[<^11 “ ^^*<^ 22 ] 


[a22 — 


where 


\ u 


E for plane stress 


for plane strain 


(2 17) 
(2 18) 
(2 19) 

(2 20 ) 
(2 21 ) 
(2 22 ) 

(2 23) 
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V for planestress 


= < (2 24 

( for planestram 

2.2 Linear elastic fracture mechanics 

Below, a brief outline of the theory of fracture mechanics for linear, plane elasticity, is giver 
Consider a vertex A, as shown in figure (2 1 ) In the case of planar elasticity, the Airy sires 



Figure 2 1 An arbitrary domain with a re entrant corner The local cartesian and polar coordi- 
nate system with respect to the A^ertex A 


function U satisfies the biharmonic equation 

/ 5^ 15 1 5^ \ / 5^ 15 ^ \tj _ ^ 

\ 5r^ r 5r \ 5r^ r 5r 50^ / 

Let us consider solutions of the form (see [ 12 ]) 

U = U(r,9) = r^+^Fi9) 

Therefore , F(9) must satisfy 

F"" -f 2(A2 + 1)F" + (A^ -1)F = 0 


(2 25) 


(2 26) 


(2 27) 


Avhere the prime represent differentiation with respect to 9 The general solution of equa- 
tion (2 27) for A 0 and A ±1 is 

F(9') — di cos(A — 1)9 -+■ 02 cos(A -t- 1 )^ -I- 03 sin(A — 1)9 -l- (I 4 sin(A + 1)0 (2 28) 
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The modifications necessary for the cases A = 0 and A = ±1 are ob\ious The stiess components 
are related to the stress function U as follows 


r^-'[(A + l)F(^) + F"(0)] 


1 ^ 1 d'^U 

^ r dr 

= A(A + («) 


(2 29j 
(2 30) 

rdrde ' r^de ~ ^ 

Let us assume that the edges at 0 = ±Y are stress free which implies that gq = t^q = 0 
on these edges Consequently, from equation (2 28), (2 31) and (2 31), after straightforward 
algebraic manipulation, we have 


cos(A — 

l)f 

cos(A -h 1 )y 

0 0 


f \ 

ax 

Asm(A - 

-l)f 

sin(A + l)f 

0 0 

< 

<^2 

> 

0 


0 

sm(A — 1 )y sm(A-l-l)Y 



0 


0 

AC0S(A — l)y COS(A + l)y_ 


. ^4 , 


0 


(2 32) 


where A is defined as 


A- 1 
A + 


(2 33) 


This gives rise to two sets of decoupled simultaneous equations in terms of ai, and as, 04 
A non trivial solution exists only if the corresponding determinant vanishes This is possible if 
either 


cos(A — 1)^ sin(A + 1)^ — Asm(A — 1 )^ cos(A + 1 )— = 0 (2 34) 

Z L Z L 


or 

sin(A - 1 )| cos(A + 1)| - Acos(A - 1 )| sin(A + 1)| = 0 (2 35) 

To find non trivial solution for equation (2 32) we set 03 = 04 = 0 and find A such that equa- 
tion (2 34) satisfied, or we set ai = as = 0 and find A such that equation (2 35) is satisfied 
Equations (2 34) and (2 35) can be simplified to 


sin Aa + A sin a = 0 (2 36) 

sin Aa — A sin O' = 0 (2 37) 


Now A can be complex and either a simple or multiple root of equation (2 36) or equation (2 37) 
If A is complex, then its conjugate is also a root of equation (2 36) or equation (2 37) , and both 
the real and the imaginary parts of are solutions of equation (2 25) If A is a solution, 

then -A is also a solution, however the corresponding stress field has finite strain energy only if 
the real part of A is greater than zero 
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Assume that ^ (i— 1,2 ) is a solution of equation (2 36) and Ap^ is real and simple Then 

from equation (2 32) 


Oi 


Gi cos(Ap^ - 1)| + 02 cos(Ap^ + 1)| 
-l)| + a2Sin(Al^^ + l)^ 


sm(A 


0 

0 


Let us define 


Q 


A) _ 


cos 


(a;‘> - i)f 


A[‘’sin{A™ - 1) 


( 1 ) 


2 


(2 38) 
(2 39) 

(2 40) 

(2 41) 

(2 42) 

sm(Af^ + l)f Ap^ + 1 cos(Af^ + l)f ^ ^ 

Equation (2 41) is symmetric with respect to 9 and equation (2 42) is anti symmetric with 
respect to 0 From equation (2 41) and equation (2 42), expressions for stress and strain can be 
derived The stress and displacement field corresponding to equation (2 41) are called Mode 1 
fields which are given by 


cos(aJ^^ + l)f sin(Aj^^ + l)f 

With this notation, equation (2 26) can be written as 

U = A''^+i(cos(Af ^ - 1)9 + Qf ^ cos(aS^^ + 1)9) 

Similarly, if a[^^ is a real, simple root of equation (2 37), then the stress function is 


U = r^*''+^(sm(Ap^ - 1)^ + Qf ^ sm(Af^ + 1)^) 


AO 


,( 2 ) 


A2) 


where is defined as 


( 2 ) _ sm(Ap^ - l)f _ Ap^ - 1 cos(Af^ - 1) 

Vi , — /ON /2) 


u 


U 


( 1 ) 

'll 

( 1 ) 


’[(k - Ol‘>(A,(l) + l))cosAf>9 - Af>cos(A™ 
1 _A.“>, 


2)9) 


2G 

which can be written as 


+ Q['H\{1) + 1)) sin A;'^0 + Af ) sin(A!'^ - 2)9] 


( 0 , 


AO 


AO 


u 


(0 - 


'1% 

2i 




where k is Kolosov’s constant given by 


AC = 3 — 4i/ for plane strain 
K = - — - for plane stress 


(2 44) 
(2 46) 

(2 46) 


(2 47) 
(2 48) 


1+F 

The Mode I stress tensor components are 

-w = A«rA!‘’-'[2-0«(A‘’ + l)co8(A«-l)9-(A«-l)cos(A«-3)e] (2 49) 

S = A«rA.“’-'[2 + <3<‘>(A('> + l)cos(A«-l)e+{A«-l)cos(A«-3)9] (2 50) 

SS = A.<‘VA.‘’-*[(a['> - l)sin(Af' -3)9 + (3*‘>(A<'’ + l)sm(Af> - 1)9] (2 51) 


a 
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The stress and displacement fields corresponding to equation (2 42) are called Mode II fields 
which are given by 

4? = “ ^fV^(2) + 1)) sm - Af’ sin(Af) - 2)6] (2 52) 

4? = + <5fHA,(2) + 1)) cos + Ap’ cos(Af) - 2)6] (2 53) 

which can be written as 

I " (2 54) 

The Mode II stress tensor components are 

^Jil = [2 - Qf ) (Af ^ + 1) sin(Af ^ - 1 ) 6 - (Af ^ - 1) sm(Af ^ - 3)0] (2 55) 

45 = [2 + gf ) (Af ) + 1) sm(Af ^ - 1)0 + (Af ^ - 1) sm(AP’ - 3)0] (2 56) 

- l)cos(Af) - 3)0 + gf^(Af) + 1) cos(Af) - 1)0] (2 57) 


Thus in the general case m the vicinity of the re-entrant corner, the solution can be written using 
equation (2 46) and (2 54) as 


{u}=E 



(2 58) 


The coefficients and are called generalized stress intensity factors In case of a crack we 
have Q! = 27 r Equations (2 36) and (2 37) are identical (sin 27r A = 0) , and all roots are real 
and simple 

A<'> = A« = ±i, ±1 ±2. ±1 (2 59) 

The coefficients Af ^ and Af ^ are related to Mode I and Mode II stress intensity factors of 
linear elastic fracture mechanics In linear elastic fracture mechanics three modes of crack opening 
are present in three dimension In plane elasticity, only first two modes are present These are 


1 Mode I - Tensile Mode or Opening Mode 

2 Mode II - Shear Mode or Sliding Mode 

Below, we give the expression for the displacements and stresses for the mentioned two modes 


Mode I 

For Mode I type of loading 


the solution in the vicinity of the crack tip is given by 
Ki 


an = 


a'22 — 


, :COS- 1 

\/27rr 2 


an 


6 30, 

sin - sin 
2 2 ^ 


Ki 0,, 0 30, 

^cosjll+smjsm-] 

Ki 0 0 30 

cos -sin 2 cos ^ 


(2 60) 
(2 61) 
(2 62) 


10 



where Kj is the stress intensity factor in Mode I The displacement components ui and uy are 
given as follows 


Ki(T\k 

^1 = — 1 — ) cos 


e_ 

2L 


ti2 


G \2Tr 
G V27ry 


1 — 2u + sin^ ^ 


sin - 
2 


2 — 2v cos^ 


0- 

2 . 

2 . 


(2 63) 
(2 64) 


<^33 = 0 for plane stress 

<^33 = + (722) for plane strain 


Mode II 


For mode II type of loading the solution in the vicinity of crack tip is given by 

(Til = 


5(2 + cos 5 cos-] 


0'22 


(2 65) 
(2 66 ) 

where Kn is the stress intensity factor m Mode II The displacement components Ui and U 2 
are given as follows 


2 

Kjj 9 9 39 

—p= sin - cos - cos — 

2 2 2 

Kji ^ 9 39^ 

(Ti 2 = — ;==cos-[l - sm-sm— ] 


Kn / ^ V ^ 

“”2L 


2 — 2u + cos^ - 

2J 


Kn / r \ I 9 r .9-] 

= ■^( s ) iJ 


(2 68) 

(2 69) 


As it IS evident from the functional form of that singularities m the form of are present 
This means that when r — 0, (Jn and CT 22 becomes infinite 


Mixed mode analysis 

The near tip displacement fields for combined Mode I and Mode II loading are obtained by 
superimposing the two fields and is given by 

9 r - „ O 


Ki / r \ h 9, I , n 2^1 Kn 

Ki / r \\ 9 , , „ 2^1 

= SgU =“ 2 ''= + ^“^“® 2 )" 

Similarly, the stress field is given by 


( — 
1^2 


2G + 

^11 ( T \\ 9 . 



Ki 


cosf ^1 -sin|sm^j 




Trr 


9 f 9 

cos I ( 1 + sin I sin 


ft ft ^9 
cos I sin I cos ^ 


¥) 


+ 


K 


II 


V2 


Tir 


sin f f2 + cos I cos ^ 

a ft ^ft ^ 

sm ^ cos I cos ^ 


9 


9 39' 


cos ( 1 — sm sin 2 


(2 70) 
(2 71) 


(2 72) 
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The above general definition of the asymptotic solution near the crack tip is fundamental to 
the process of quasi static crack propagation Below, we define the J integral that is used to 
characterise the criticallity of an existing crack 

2 2 1 J integral 

Taken along a material curve F which is traversed in the counter clockwise direction between the 
two crack sides , the J integral is defined as 

J = J^(Wdx2 - T^^^ds) (2 73) 

where W = o^^dsij is the deformation energy per unit volume, s is the arc length and % is 
the traction exerted on the body bounded by the curve F and crack surface 

T) — O' Tlj 

where Uj is the unit outward normal to the curve F Two important characteristics of J integral 


^2 



Figure 2 2 Computation of J integral 


are 


• It IS path independent for linearly elastic bodies 

• For linearly elastic bodies, J integral represents the energy release rate and is same as the 
energy release rate G, proposed by Griffith [15] 
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It can be shown for Mode I loading condition 


and for Mode II loading condition 


Thus for mixed mode analysis, 


.=3 

E* 


1-!^ 

^ E* 


j = M + iSt 
£• £• 


(2 74) 


Thus we note that the computation of the stress intensity factors is essential for obtaining the J 
integral 


Computation of stress intensity factors 

The stress intensity factor can be computed by several methods Below we consider one of 
the method outlined m [3], which employ contour form of the interaction integral The local 
coordinates are the crack tip coordinates with the Xi axis parallel to the crack faces (see figure 
2 2) Two states of cracked body are considered following [3] State 1 corresponds 

to the actual state and state 2 is an auxiliary state which will be chosen as the 

asymptotic field for Mode I and Mode 11 The J integral for the sum of the two states is 


j(l+2) _ f _(^1^ ^ ^ 


rijdr 


Expanding and rearranging the terms give 


ja+ 2 )^ J(l) + j( 2 )+j( 1 . 2 ) 

where is the interaction integral for states 1 and 2, 


ja.2)= f [p^(i. 2)5 ^^a)M_cr;2)^ln,dr 

JpL dxi dxii 


and TT(^>2) is the interaction energy, given by 

TT(b2) ^ ^ 

Substituting equation (2 74) for the combined states and rearranging the term gives 

j(i+2) = J(1) + J(2) + + Kf^Kf}) 

Equating equation (2 76) and (2 79) leads to the following eqaution 


(2 75) 


(2 76) 


(2 77) 


(2 78) 


(2 79) 




(2 80) 
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Now choosing state 2 as Mode I asymptotic field with Kf'^ = 1 and Kf^ = 0 gives Mode I 
stress intensity factor for state 1 or actual state in terms of interaction integral as 

= KjUMode /) ^2 81 ) 

2 

Then choosing the auxiliary state as that for Mode II asymptotic field with = 0 and 

(2) 

Ku — I gives Mode II stress intensity factor for state 1 in terms of interaction integral as 

= EljihMode II} ^2 82 ) 

2 2 2 Crack growth criteria 

We briefly review the crack growth law used to specify the direction of crack growth The criteria 
for determining the crack growth direction can be summed up m the following three categories 

1 The maximum energy release rate criterion 

2 The maximum tangential stress criterion 

3 The minimum strain energy criterion 


In the present work, we used the maximum tangential stress criterion or the maximum circum- 
ferential (hoop) stress criterion, which is identical to maximum energy release rate criterion for 
linear elastic fracture mechanics The maximum tangential stress criterion states that the crack 
will propagate from its tip in a direction 0c so that the circumferential stress agg is maximum 
Therefore, the critical angle dc deflnmg the radial direction of propagation can be determined by 
setting the shear stresses associated to zero because this corresponds to the principal direction 
Thus 


= cos - -Kf sin 9 + -Ku (3 cos 0 - 1) = 0 (2 83) 

27rr 2 L2 2 J 

This leads to the equation deflnmg the angle of crack propagation (in the tip coordinate system) 


do 


Ki sin 9c + Kiiifi cos 0c — 1) = 0 


(2 84) 


Solving the above equation gives 

e, = 2arctani(;g±y(£7+8) 


(2 85) 


For Ku = 0, 0c = 0 , 1 e the crack propagates straight, without turning 
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Chapter 3 

Finite element formulation 


3.1 Finite element formulation using energy principle 

The total potential for the plate is given by 

N 

n (u) = (u), (3 1) 


e=l 


where tt® is the total potential of the non-mtersecting (but adjacent) sub-domains e which are 
part of the domain (N sub-domains are considered here) The total potential can be expressed 
in terms of internal strain energy {7^®^ and external work done as follows 

TT® (u) = 17^®) - TF(®) (3 2) 

n(u) = ^ / {(^xx{^)£xx{n)+crxy{n)'y^y{u)+ayy{n)eyy{v.)}dA- jj^fxU+fyv)dA- J^{TxU+Tyv)ds 

(3 3) 

where {f } is the body force vector, {T} is the traction vector and {u} is the displacement vector 
given by 


= { A } 
= {?:} 

Minimizing the total potential energy gives 


d^^)(n(u)) = hma-^o 


n(u + aw) — n(u) _ 


= 0 


a 


where {w} is the perturbation vector 


{w} 


Wx 

W,, 
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5^^^(n(u)) - Jj,<7xxin)ea:x{w) + a-a;j,{u)7(w) + <Jyyin)eyy{w)) - J {fxWx + fyWy)dA - 


{Tx'^X TyVJy'^dS 


Hence the bilinear form can be written 


^(u,w) 


(0-a:i(u)£ix(w) + axy{vL)jxyi'^) + (^yy{^)£yy{^))dA 


-P’('W) = J {fxU!x + fyWy)dA + j^(TxWx +TyWy)dS 


3.2 Finite element approximation 


{ufsm} = 


EILi 


{ufsm} = 


$1 0 $2 0 *^3 0 

0 $1 0 $2 0 $3 


= [^]{U} 


Similarly, 


f Ei=i i 

{w} = 

I Er=i%^^^ J 

{w} = [#]{W} 

Triangular elements are used in the finite element approximation employed m this study, along 
with hierarchic shape functions of order p (p < 3) n is the number of degrees of freedom for 
each element and it is equal to for triangular elements are shape functions in global 

coordinate system x and y 

Let the stress and strain vectors corresponding to {u} be {a} and {e:} which can be defined as 


{a} = < cr 


Prom the generalized Hooke’s law which relates stress components to the respective strain com- 
ponents in global coordinate system, we have 

W = [Q]{e} (3 7) 
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wh6r6 ths rnatsnal stiffnsss matrix [Q] is given by 


[Q] = 


E vE n 

l-i/2 1-1,2 U 


uE E 


0 


1 — 1/2 1 — 1/2 

0 0 G 


for plane stress 




[Q) = 


i^E 


1— i/— 2i/2 1— 


uE 


E(l-iy) 


0 


0 


1 — z/— 2i/2 1— i/— 2i/2 

0 0 G 


for plane strain 


The strain can be represented as 

{e} = [D]{u} 

where [D] is a differential operator in terms of global coordinates such that 


[D] = 


0 


0 f 

ay 


A A 

L dy dx J 


Hence, the strain field using finite element approximation is given by 

= [D]{ujir£;A^} 

Substituting the expression for we have 

= [B]{U} 

where [B] is given as 


a<l>i 

0 

d^2 

0 

a^?. 

0 

- 

dx 

dx 

dx 


0 

dy 

0 

a^2 

ay 

0 

dy 









- dy 

dx 

dy 

dx 

dy 

dx 

- 


[B] = 


Similarly, for 

£(w) = [B]{W} 

Thus, in finite element approximation the bilinear form can be written as 

Pi^FEM,^) = / {e'^i'^)}{cr{uFEM)}dA 

Ja 

/3(u™„,w) = [B]{W})^[D]([B1{U} 


(3 8 ) 


(3 9 ) 


( 310 ) 


( 311 ) 


( 312 ) 


(3 13 ) 


(3 14 ) 
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P{nFEM,w) = {W}^[B]^[D][B]{U} 
The linear functional can be written as 

F('w) = f fwdA + f Twds 


= Jj.^f{^}dA + J^Jwy{T}ds 


(3 15) 


F(w) = {W}^ / + / [$]'^{T}ds 

^Ja Jt ■* 

Now, substituting the above expressions for bilinear and linear functional , we get 

/3(uj?em,w) = F{-w) 

{Wr[^[B]'"[D][B]d^]{U} = {W}^[Jj_^f{f}dA + J^[^f{T}ds 

Since, {W} IS arbitrary, we have 

[K]{U} = {F} 

where [K] is the element stiffness matrix given by 


[K] = / [B]’’(D][B]d^ 
J A 


and {F} is the load vector given by 


{F}= f [$f{f}dyl+ f [^f{T}ds 
J A JT 


(3 16) 


(3 17) 


(3 18) 


(3 19) 


3 2 1 Computation of stiffness matrix and load vector 

After expanding the terms of equation (3 18), the entries for element stiffness matrix can be 
written in the following form for a given i and j 


K,„. 


(2.-l)(2,-l) 


(On* 3,x^i,x + Q33^j^y^i^y)dA 


K(2 i- 1){23) = / {Ql2^j,y^t,x + j,x^i,y)dA 

Ja 

if'(2t){2j-l) = L iQl2^],x^i,y + Q3Z^j,y^i x)dA 


^(2t)(2j) 


{Q22^o,y^i‘,y ■f' Q^^^o,x^i,x)dA 


Here % and j vanes from 1 to NDOF (where NDOF is the number of degrees of freedom) 
Similarly, from equation (3 19) we get 


F( 2 i-i) = J fx^idA + J^Tx^tds 

F{2i) = J fy^idA + Ty^ids 
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3.3 Geometry mapping 

Numerical integration of the stiffness matrix entries and the load vector entries necessitates the 
transformation of the physical coordinates (x, y) to the master element coordinates rj) Hence 
instead of defining the shape functions in the physical coordinate s} stem, it is defined in master 
element coordinates denoted by The geometry is expressed in terms of shape functions 

3 

(3 20) 

2=1 
3 

(3 21) 

2=1 

where and y® are the physical coordinates of the i-th node of the element e and is the shape 
functions defined m master element coordinates (i^, y) 

Linear mapping 





Figure 3 1 Linear mapping 

Linear mapping is generally used when the sides of the physical elements are straight edges 
Then the shape functions are given as 

Ni = l-(-V 

iV2 = e 

Ns = 77 

Substituting the above relations, we get 

X = - C - T]) + xl^ + xlr) 

y = yt{i - ^ - v) + y2^ + ylv 

or 
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f X — xf 
{y-Vi 


{xi - xi) {xi - xf) 1 r ^ 1 

(yl - yt) (yl - yf) J 1 y j 


V 


1_\ (y| - y?) 

|J| L -(y|-yO 


-(x| - xf) 
(x| - xf) 



where [J] is the jacobian matrix given by 


[J] = 


(x|-xf) {xl-x\) 

(yl - yi) (yf - yf) 


The derivatives of the shape functions are calculated using chain rule 


^ 

dx dx dr] dx 

^ dN^dr] 

dy d^ dy dr] dy 


The derivatives for linear transformation are given by 


dx 

dy 

dr] 

dx 

dr] 

dy 


- ^5) 


(3 22) 


(3 23) 


(3 24) 
(3 25) 
(3 26) 
(3 27) 


Numerical integration 

The numerical integration over a general triangle is achieved by mapping the triangle to a defined 
master element and then performing the area integral The mapping can be either linear or 
quadratic or any order which depends on the geometry of the general element is described in the 
above paragraph After transformation, the integral value is given by 

[ F{x,y)dxdy = f F{^,r])\3\d^dr] 

Ja^ Ja 

. MINT 

2=1 

where A® is the area of the physical element, A is the area of the master element and NINT is 
the number of integration points 


20 



Chapter 4 

Generalized finite element formulation 


4.1 Introduction 

In the standard finite element method, the approximate solution is constructed by employing 
piecewise polynomial function of degree p A major drawback of this method is that for the 
domain with re-entrant corner (or crack), severe mesh refinement, near the vertex of the corner 
is required Thus the size of computational problem increases This is because the solution is 
very unsmooth m the vicinity of the crack tip However in case of linearised elasticity based 
fracture mechanics, the representation of the unsmooth part of the solution is known m the 
vicinity of the crack tip, as a series Hence, incorporating this special function in the finite 
element representation will significantly improve the quality of the finite element solution in 
the vicinity of the crack tip Thus the necessity of the mesh refinement would be drastically 
reduced for this type of finite element analysis Further, the critical parameter (for example the 
J integral) can be obtained more accuartely as compared to conventional finite element analysis 
Below, a detailed outline of the generalized finite element method (GFEM) is given 

4.2 Construction of the conforming generalized finite el- 
ement approximation 

In this section, we describe how the finite element approximation can be generalized to include 
special functions The standard shape functions used in finite elements result in mapped poh- 
nomials and are known to approximate smooth functions well However, if the solution is know n 
to be harmonic, it would be advantageous to use harmonic polynomials to form the basis But 
the problem of continuity of approximation or conformity should be addressed properly to make 
sure that the approximation has finite energy Conformity has been addressed in the standard 
finite element method (FEM) by constructing element shape functions in such a way that the 
proper number of degrees of freedom match at the boundaries of neighbouring elements and thus 
create a continuous approximation Using any of the standard FEM elements and an associated 
family of shape functions of degree p, we can exactly represent all polynomials of degree p These 
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shape functions are pasted together from neighbouring elements and form a basis by suppressing 
any linear dependencies, e g , rigid body modes in elasticity The resulting system of equations 
are solvable The partition of unity method allows creation of conforming approximation using 



Figure 4 1 Linear support hat function 



Figure 4 2 Special function to be used for enrichment for e g , ® 



Figure 4 3 Generalized basis function given by ®(1 - 2x) over the domain (0, 0 5) 

special functions Given a set of overlapping subdomains or patches {GJ (see figure 4 1) and a 
set of functions with desirable approximation properties associated with each patch = {ip]}, 
(see figure 4 2) we define the PUM approximate solution over the domain Q. (see figure 4 3) as 

upt/M = (^]'^]) (4 1) 

^ j 
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where ‘ip^^{x,y) is a special function from the set {ip^} The space Tj is called the patch space 
The are the constant coefficients The sequence of functions {#j} (one for each patch of the 
set IS a partition of unity on Q and serves to enforce contmuitj- 

Although the partition can be constructed in any manner as long as it satisfies the definition, 
in the present work we construct it using linear hat functions over coarse patches of triangular 
elements (see for one dimension example figure 4 1-43) The patches and each associated 
member of the sequence {$j}must satisfy the conditions stated belo'v\ 


{Oj IS 

an ( 

3 pen co\er of Q, } 



e 

0 

0 

< 

(4 2) 

sup$i 

G 

dosure{Q.i) V i, 

(4 3) 

E*. 

= 

1 on f 2 . 

(4 4) 

114.11? 

< 

Coo, 

(4 5) 


The displacement field, at a crack tip, is contained in the span of the following four functions 


( 9 6 9 9 ') 

{ 7 (r, 0)} = |\/f cos -, y/rsm -, Vr sm - sm9, ^/rcos - sin^j 


(4 6 ) 


These functions will be used to enrich the standard finite element shape functions to give better 
approximation Although more general form of enrichment function can be used, we used the 
above span to enrich the shape function used by Belytschko et al [4] Thus we can completely 
generalize the finite element approximation, allowing any family of special functions to be mixed 
with the mapped polynomial shape functions The generalized finite element approximation is 
given by 

rivert npEM 

^GFEM = ^ ^. 7 ^.?) + '^kNk (4 7 ) 

z=l j k 

where rivert is the number of hat functions at element vertices, is the number of special functions 
in the patch space associated with the vertex and ufem is the number of finite element side 
modes and interior modes 

Note that the traditional FEM basis functions are evaluated on the master element (^, 77 ) but the 
higher order special functions are expressed in terms of unmapped physical coordinate system of 
the problem (x, y) 


4.3 Linear dependency 

As mentioned earlier, the GFEM may lead to a singular stiffness matrix because of linear depen- 
dence of the local approximation functions On the contrary, in traditional FEM the solution 
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corresponding to zero eigen value is removed by imposing essential boundary conditions However 
this IS not possible m GFEM because the eigen functions associated with the zero eigen value are 
not known generally To solve such a linear system of equations one can use several methods like 
the perturbation approach [6] In the present study, perturbation technique is incorporated in 
the solver code to solve such system of equations, though singularity does not arise in our case 
When [A] is a singular matrix, then a particular solution can be obtained by 

[A]e = [A] + e[I] 

where s is the perturbation, usually taken very small and [I] is the identity matrix Thus [A^- 
is positive definite matrix and can be inverted by any standard method 


4.4 Adaptive integration 


Since in the present work we are enriching the usual finite element space with local asymptotic 
solution in the vicinity of crack tip which have singular derivatives at r = 0, we must take care of 
numerically integrating their entries in element stiffness matrix or load vector This is done using 
adaptive gauss quadrature integration, in the elements where the special functions are employed 
For all other elements the standard gauss quadrature rule is applied In the following paragraph, 
adaptive integration scheme will be discussed elaborately 
The GFEM approximation will be defined as 

1 9 

M^,v) cos -Nm{(,v) 

where v) the traditional finite element shape function defined in the master element and 
IS the new auxiliary shape function Now to compute the stiffness entries, we need to 
find the derivatives of the auxiliary shape functions which we can find by applying chain rule 


dm,rj) ^ dNm{^,v) d^ dNmi^,ri) dr} 
dx l dx dr] dx 


7-2 cos 2 + ^rn 


. d , 1 6 

(C,^)^(’'^cos-) 


So, the stiffness entries corresponding to the auxiliary degree of freedom involves singular deriva- 
tives To compute the stiffness entries for such row and column which involves the auxiliary 
degree of freedom, we divide the master element into four equal parts and compute the area 
integral on each of the sub master element and added to give the total value for that element 
(see figure 4 4) This value is then checked with the value calculated in the previous iteration, if 
the error is within the range of the tolerance specified in the code then no further subdivisions 
are carried out and calculation of next stiffness entry is carried out However, if the error is more 
than the tolerance specified then out of the four sub master element, the sub master element 
asociated with the crack tip is redivided into four equal triangles and the integrand is evaluated 
on each of the element parts and summed to give the total value (see figure 4 4) Special care 
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should be taken to ensure that the weighting factor of the standard gauss quadrature rule gets 
modified after a triangle is divided For example, m the first iteration after division of the master 
element into four equal triangles the weighting function for all the sub master elements should 
be divided by four While the situation in any other iteration greater than one is more complex, 
since it involves both the sub master elements obtained m the previous iteration and the integral 
domain obtained after the division of the triangle associated with the crack tip into four more 
parts Hence special care should be taken so that all the previously divided elements which have 
not been subdivided should retain their weighting function pertaining to the previous iteration 
and for the triangle which has been sub divided, its weighting function should be divided by four 
again to get the weighting function for each of the sub divided triangles In the present work a 


3 Ennched Node 



Figure 4 4 Scheme for Adaptive Integration 

ninth order integration rule is applied which results m nineteen integration points In each of the 
sub divided triangles, the same order integration rule is applied to maintain the simplicity of the 
code An important aspect of the present adaptive integration procedure is that we apply the 
prior knowledge of the node corresponding to crack tip which has been enriched and subsequently 
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attack the node to achieve desired accuracy This approach ensures that the stiffness entries are 
accurate However one has to be cautious since the same integration rule is applied to each of the 
sub divided domain, proper care should be taken to transform the sub master element coordinate 
to the actual master element This is due to the fact that the shape functions are defined in 
actual master element and not in the sub master element This is taken care of by a linear map 
of the sub master element coordinates to the actual master element coordinates since the master 
element is an isoceles right angled traiangle having straight edges 

An important observation is that the present approach takes longer time for diagonal elements 
This IS due to the fact that the off-diagonal elements exhibit weaker singularity than the diagonal 
elements In this work, we imposed a tolerance of Note that, we could have carried out a 

single computation of unsmooth integrand, and obtain the desired level of subdivision required 
in the master element, for elements at the crack tip This can then be used to do the integration 
in one shot, for all the elements at the crack tip 

4.5 Solution technique 

The solution phase of the finite element computer programs consists of three parts 

1 Assembly of the stiffness matrix and load vector from the element level stiffness matrices 
and load vectors, 

2 Enforcement of the principal boundary conditions, and 

3 Solution of the system of simultaneous equation 

In engineering problems the size of the system of equations is generally quite large, and the entire 
global stiffness matrix cannot be stored m runtime memory, even if only the non-zero entries are 
stored It is important to design the solution algorithm so that the arithmetic operations in the 
solution phase can be performed as efficiently as possible 

In the present work, we have used the frontal solution method The frontal solution method was 
proposed by Irons [12] Its mam feature is that the order m which the gaussian elemmation is 
performed is element oriented rather than node oriented It is not necessary to number nodes for 
minimum bandwidth Another important feature of this method is that assembly of the stiffness 
matrix and the gaussian elimination is carried out concurrently 
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Chapter 5 

Numerical results for stationary crack 


In the present chapter we will discuss the results obtained Both straight and arbitrar\ shaped 
crack segments were considered Initially a comparison is drawn between FEM solution and 
GFEM solution to prove the efficiency of the proposed methodology The cases w^ere solved using 
hierarchic shape functions of order p (p<3) This could be very easily extended to even higher 
order approximating polynomial In all the cases, the enrichment functions used were conforming 
to the one proposed by Fleming et al [4] However more general enrichment functions could be 
employed 

Crack modelling 



Actual Crack 



Figure 5 1 Crack modelling 

In the present work, any arbitrary shaped crack is modelled as a series of straight line segment 
connected at nodes as shown in figure 5 1 To ensure the discontinuity along the crack faces, two 
nodes were created at the same position one on the upper surface of the crack and the another 
one on the lower surface This ensures that the element discretization above and below the crack 
faces are independent of each other The local coordinates at the crack tip are taken parallel to 
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the end line segment (see figure 5 1) 

In all the cases the full plate analysis is done to elaborate the computational strategy in case of 

multiple singularities in the domain For all the problems, plane strain condition is chosen and 

the following material properties, boundary conditions and loads are taken 

Material properties E = 100 kvsi and = 0 3 

Boundary conditions Four corner nodes fixed 

Traction = 1 0 psi on the opposite edges 




Figure 5 2 Boundary conditions and loading conditions 
The plate is taken to be of dimensions 5m x 3 5m (as shown in figure 5 2) 

Problem 1 Centrally located straight edge small crack 

A center cracked plate with the following parameters is analysed 

Domain = 5 x 3 5 
Crack Length = 0 5 in 

Two different mesh is employed, one is coarser(see figure 5 3) and the other one is finer(see figure 
5 4) The computed values of the various crack parameters using either the conventional FEM 
or GFEM are reported in table 5 1 and 5 2 Note that in the computation of stress intensity 
factors, the contour is taken to correspond to the outer boundaries of either one layer of elements 
adjacent to the crack tip or two layer of elements 
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Figure 5 3 Straight edge crack of length 0 5in with coarse mesh 



Figure 5 4 Straight edge crack of length 0 5in with fine mesh 
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Table 5 1 Stress intensity factors computed for coarse mesh 


Scheme 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 
P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Kj m psiVin 

Kjj in psiVm 

FEM 

634 

1 

1 

0 871998271 

-0 00186705222 

0 871918492 

-0 00213316604 

GFEM 

634 

1 

1 

0 764668694 

-0 0017483637 

0 766995734 

-0 00157027579 

FEM 

634 

2 

1 

0 82840364 

-0 000819087014 

0 828413002 

-0 000812332551 

GFEM 

634 

2 

1 

0 978282491 

-0 00123810708 

0 978622215 

-0 00106178627 

Tabl 

e 5 2 Stress intensity factors computed for refined mesh 

Scheme 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Kj: in psi\/in 

Kj} psi\/in 

FEM 

2564 

1 

1 

0 626637336 

-0 000337844132 

0 626767391 

-0 000283837087 

GFEM 

2564 

1 

1 

0 500493715 

-0 000460228968 

0 50052228 

-0 000217512687 

FEM 

2564 

1 

2 

0 849307746 

-0 000258345046 

0 849500731 

-0 000256760112 

GFEM 

2564 

1 

2 

0 855961206 

-0 00104771529 

0 855961206 

-0 000787968813 


Infinite Plate in Mode I Loading 


Ki = 0 895 


psiv m 


Kji — 0 psiyin 


Discussion 

1 The theoretical value of Ki and Kn for finite plate, after applying finite plate corrections, 
are Kj — 0 895 and Kn = 0 The value obtained by GFEM and FEM with either 
p = 2 and one layer or p = 1 and two layer are sufficiently close to the theoretical 

values 
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2 The contour integral method applied in the study is highly sensitive to the distance of the 
contour from the crack tip As the contour is brought closer to the crack tip, the quaht} 
of the extracted stress intensity factor deteriorates 

3 The stiffness intensity factors obtained at the left crack tip and right crack tip do not match 
exactly because of the lack of symmetry m the mesh used Though, it should be noted 
that the discrepancy is insignificant 

4 For coarse mesh, the results are not good beacause the crack length is too small The 
contour integral employed in the extraction of the stress intensity factors will be severely 
affected by the neighbouring crack tip This is the reason that Fleming et al [4] suggested 
mesh refinement at the crack root for getting better results 

Problem 2 Centrally located straight edge long crack 

In this problem a centre cracked plate with the following parameters is analysed 

Domain = 5x35 
Crack Length = 10 in 

A. different crack length is chosen to demonstrate the effect of increasing crack length on the 
solution and how the higher order shape function affects the solution In this problem the mesh 
employed is coarsh(see figure 5 5) The computed values of the stress intensity factors are reported 
m table 5 3 and 5 4 



Figure 5 5 Straight edge crack of length of lin 
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Table 5 3 Stress intensity factors computed for p = 1 


Scheme 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Kj m psi\/m 

Kij psi^/m 

FEM 

628 

1 

1 

0 90438288 

-0 00278504271 

0 904393975 

-0 0028196911 

GFEM 

628 

1 

1 

0 744402238 

-0 00184098385 

0 74466355 

-0 00191254745 

FEM 

628 

1 

2 

0 764668694 

-0 00236289666 

0 766995734 

-0 0023754772 

GFEM 

628 

1 

2 

1 25931867 

-0 00262408054 

1 26092989 

-0 00329535918 


Table 5 4 Stress intensity factors computed for p == 2 


Scheme 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Kj in psi\/m 

Kii psiy/in 

FEM 

628 

2 

1 

0 841305802 

-0 0017483637 

0 841326947 

-0 00157027579 

GFEM 

628 

2 

1 

0 969028848 

-0 001050745056 

0 968735869 

-0 00149841222 

FEM 

628 

2 

2 

1 21598151 

-0 0017483637 

1 21636486 

-0 00157027579 

GFEM 

628 

2 

2 

1 27639092 

-0 00104771529 

1 27570483 

-0 000787968813 


Infinite plate m Mode I loading 


Ki = 12B> ps%y/%n 


Kii = 0 pse V in 


Discussion 

1 For p = 1, when the contour corresponds to two layer of elements, GFEM is very accurate 
while the standard FEM is not Further using one layer of elements around the crack tip 
seems to be inefficient for the extraction of stress intensity factors 

2 For p = 2, GFEM is clearly far superior to standard FEM 


32 























3 Note that in this problem, the interference from the neighbouring crack tip is very small 
leading to more accurate extracted quantities 

4 Inspite of a coarse unsymmetric mesh, the left and right stress intensity factor have in- 
significant mismatch 


Problem 3 Arbitrary shaped crack 


A.n arbitrary shaped crack with the following parameters is analysed 

Domain = 5 x 3 5 
Crack Segment 

nodel = (2 00, 1 70) 
node2 = (2 25, 1 70) 
nodes = (2 50, 1 80) 
node4 = (2 75,1 80) 


The plate with the crack is shown in figure 5 6 The computed stress intensity factors are reported 
in table 5 5 



Figure 5 6 Problem 3 


Table 5 5 Stress intensity factors computed for different p 


Scheme 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Kj in 'ps%\fm 

Ku psiy/in 

GFEM 

636 

1 

1 

0 990401579 

0 0154856607 

0 982036791 

0 010429162 

GFEM 

636 

2 

1 

1 2372214 

0 0320909648 

1 24435412 

0 0309407329 

GFEM 

636 

3 

1 

1 20827447 

0 0536397577 

1 19315666 

0 022451595 


Discussion 


1 The stress intensity factors calculated for p = 2 and p = 3 are quite close and seems 
to converge (i e discrepancy is less than 5 pc) 

2 Note that the Kn is essentially zero as compared to Kj and the crack will tend to grow 
straight The left and the right stress intensity factor do not match because mesh is crude 
and unsymmetric, however, the discrepancy is minimal 

3 To improve the result, one or two level of mesh refinemnet near the crack tip is desirable 


Problem 4 Arbitrary shaped crack 

An arbitrary shaped crack with the following parameters is analysed 


Domain = 

Crack Segment 

5 X 3 5 zn' 

nodel = 

(2 00, 1 70) 

node2 = 

(2 25, 1 72) 

node3 = 

(2 50, 1 75) 

node4 = 

(2 75, 1 70) 

node5 = 

(3 00, 1 75) 

The arbitrary shaped crack is modelled as shown m figure 5 7 1 
factors are reported m table 5 6 and 5 7 
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Figure 5 7 Problem 4 


Tal 

ole 5 6 Stress intensity factors computed for one layer 

Scheme 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Elements 
to compute 
Interaction 
Integral 

K/ in psiV^’U 

Kjj psi-s/m 

GFEM 

584 

1 

1 

1 13418414 

-0 00521811786 

1 08020361 

-0 0412967273 

GFEM 

584 

2 

1 

1 28536871 

0 0419222697 

1 28719329 

-0 063198877 

GFEM 

584 

3 

1 

1 47696526 

-0 0364505828 

1 41330589 

0 027973453 


Table 5 7 Stress intensity factors computed for two layer 


Scheme 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Elements 
to compute 
Interaction 
Integral 

Ki in ps%\/m 

Kii psi'/in 

GFEM 

584 

1 

2 

1 5457882 

-0 051235026 

1 56144341 

0 0193207251 



2 

2 

1 4236723 

-0 0542374126 

1 38937786 

0 0665384271 
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Discussion 


1 The results for 1 layer of elements is quite poor However for p = 3 the results seem to 
be close to that obtained from extraction using two layer of elements around the crack tip 

2 Again in this case Ku is negligible because the crack behaves like a centrally placed crack 


Problem 5 Crack growth for two symmetrically located crack 

The growth of two symmetrically placed straight cracks with the following initial parameters is 
analysed here and the corresponding position of the cracks are shown m figure (5 8-5 11) 


Domain 

Crack Profile 1 
no del 
node2 
nodes 
node4 
node5 

Crack Profile 2 
nodel 
node2 
nodes 
node4 
node5 


5x35 in^ 

(1 500, 1 75) 

(1 625, 1 75) 

(1 750, 1 75) 

(1 875, 1 75) 

(2 000, 1 75) 

(3 000, 1 75) 

(3 125, 1 75) 

(3 250, 1 75) 

(3 375, 1 75) 

(3 500, 1 75) 


The purpose of the present problem is to estimate the robustness of the code for moving geom- 
etry Here m each step, the crack length gets increased and the analysis is repeated The crack 
extension direction is determined according to maximum tangential stress criterion detailed ear- 
lier in chapter 2 The extracted stress intensity factors for different p values are reported in table 
58- 5 13 



step 1 



Figure 5 8 Problem 5 crack growth step 1 


Table 5 8 S1 

:ress intensity factors computed at the end of stej 

0 1 with p = 1 

Crack 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Elements 
to compute 
Interaction 
Integral 

Ki in psiVm 

Kji psi\/m 

1 

2564 

1 

2 

0 827900965 

-0 000352688633 

0 85007153 

-0 000317667391 

2 


1 

2 

0 850005244 

-0 000285569277 

0 827116411 

-0 0000478630334 


Table 5 9 Si 

tress intensity factors computed at the end of ste] 

1 with p = 2 

Crack 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Elements 
to compute 
Interaction 
Integral 

Ki m pszVm 

Kii pstVm 


2562 

2 

2 

0 853412354 

-0 000240084291 

0 869690615 

-0 000255675898 


2562 

2 

2 

0 869697323 

-0 000280393539 

2 

0 853348802 

-0 0000264159494 
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step 2 



Figure 5 9 Problem 5 crack growth step 2 


Table 5 10 Stress intensity factors computed at the end of step 2 with p = 1 


Crack 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Kj in ‘psi's/vn 

Ku psi-slm 

1 

2564 

1 

2 

0 820305226 

-0 000301340611 

0 86535299 

-0 000294438658 

2 

2564 

1 

2 

0 865398794 

-0 000570112828 

0 819049369 

0 00263806512 


Table 5 11 Stress intensity factors computed at the end of step 2 with p = 2 


Crack 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Elements 
to compute 
Interaction 
Integral 

Ki in pstVm 

Kii psiVm 

1 

2560 

2 

2 

0 848856914 

-0 000126843033 

0 883541027 

-0 000481974206 

2 

2560 

2 

2 

0 881666649 

-0 000156524428 

0 85113107 

-0 000796571682 
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Step 3 





Figure 5 10 Problem 5 crack growth step 3 


Table 5 12 Stress intensity factors computed at the end of step 3 with p = 1 


Crack 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Ki in psi^/m 

Kji psi\/m 

1 

2510 

1 

2 

0 912694996 

-0 000292969031 

0 964575631 

-0 00703674212 

2 

2510 

1 

2 

0 918500545 

-0 00261285263 

0 829219297 

-0 00526810027 


Table 5 13 Stress intensity factors computed at the end of step 3 with p = 2 


Crack 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

— 

Layer of 
Elements 
to compute 
Interaction 
Integral 

Ki m psi\/m 

Kji psi'Jm 

1 — t 


2 

2 

0 911692321 

-0 000531385793 

1 00692173 

-0 00147615579 


2566 

2 

2 

1 00926826 

-0 000218973043 

0 921226773 

0 00378581801 
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Final step 



Figure 5 11 Problem 5 crack growth step 4 


Discussion 

1 Following the theoretical analysis, geometry and location of the crack, both the crack should 
propagate along a straight line This is what we get out of our numerical analysis 

2 As the crack grows, the interaction between two cracks become more severe leading to an 
increase in stress intensity factor 

3 The crack propagation is stopped prior to merger because of the limitations of the current 
meshing process 

Problem 6 crack growth for two arbitrarily located crack 

In this problem two cracks are considered However the cracks are located off center as shown 
m the figure(5 12) The details of the domain and the crack geometry are as follows 

Domain =5 x 3 5 m^ 

Crack Profile 1 
nodel = (1 500, 1 25) 
node2 = (1 625, 1 25) 
nodes = (1 750, 1 25) 
node4 = (1 875, 1 25) 
nodes = (2 000, 1 25) 

Crack Profile 2 

40 / 

r 


nodel = (3 000 2 25) 
node2 = (3 125, 2 25) 
node3 = (3 250, 2 25) 
node4 = (3 375, 2 25) 
node5 = (3 500, 2 25) 

Each of the crack segment is taken as a series of line segments joining the above 5 nodes After 
each step of iteration, the crack profiles are given proper crack extension in the right direction 
and the whole analysis is repeated The computed stress intensity factors at the end of each step 
are reported in the table 5 14 - 5 23 The corresponding crack geometries at the beginning of 
each step are shown in figure (5 12 - 5 21) 

Step 1 



Figure 5 12 Problem 5 crack growth step 1 


Table 5 14 Stress ini 

tensity factors computed at the end of step 1 

Crack 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Kj in psi-\/m 

Kn psiV in 

1 

2584 

1 

2 

0 835341791 

0 0121879137 

0 850736133 

0 0113071723 

2 

2584 

1 

2 

0 850696036 

0 0113071723 

0 834182778 

0 0121700265 
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step 2 



Figure 5 13 Problem 6 crack growth step 2 


Table 5 15 Stress intensity factors computed at the end of step 2 


Crack 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Kr in pszVm 

Kii psi'/m 

1 

2596 

1 

2 

0 830585879 

0 00567496787 

0 851276091 

-0 00508734123 


2596 

t 

1 

2 

0 864274472 

0 00791880615 

0 828841718 

0 00307971677 
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step 3 



Figure 5 14 Problem 6 crack growth step 3 


Table 5 16 Stress intensity factors computed at the end of step 3 


Crack 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Kj in psiy/m 

Kii psi^/m 

1 

2578 

1 

2 

1 00135083 

-0 00158095561 

1 09011459 

-0 0132891375 

2 


1 

2 

1 093026525 

0 00643484579 

1 02283768 

0 00345989243 
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step 4 
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Figure 5 15 Problem 6 crack growth step 4 


Table 5 17 Stress in1 

tensity factors computed at the end of step 4 

Crack 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Ki m psiy/m 

Kij psi'/m 

1 

2570 

1 

2 

1 36818268 

-0 0211871321 

1 48479463 

-0 0986211241 

2 

2570 

1 

2 

1 4507246 

-0 0194111176 

1 36069751 

-0 000822044013 
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Figure 5 16 Problem 6 crack growth step 5 


Table 5 18 Stress intensity factors computed at the end of step 5 


Jo of 

Ele- 

nents 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Ki in psz\/m 

Kii psi-Jm 




1 71884708 

-0 0403886554 

i70 

1 

2 

1 86253646 

-0 0604739141 




1 84575045 

-0 200539113 

i70 

1 

2 

1 66395861 

-0 117210394 
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Figure 5 18 Problem 6 crack growth step 7 


Table 5 20 Stress intensity factors computed at the end of step 7 


Crack 


No of 
Ele- 
ments 


2580 


Order of 
Shape 
Functions , 
P 


Layer of 
Element to 
compute 
Interaction 
Integral 


Ki in psiy/m 


2 48713937 


2 56750505 
2 43620877 


2 44884128 


Kii pstVm 


-0 234273036 


-0 308852834 
-0 309309869 


-0 23661355 


80 


1 


2 


Figure 5 19 Problem 6 crack growth step 8 


Table 5 21 Stress intensity factors computed at the end of step 8 









Crack 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Kj m 'psi'Jm 

2 94236899 

Kn psis/m 

-0 324194315 

1 

2592 

1 

2 

2 6786357 

2 72784844 

-0 315583614 
-0 491364253 

2 

2592 

1 

2 

2 90369245 

-0 276572426 



step 9 



Figure 5 20 Problem 6 crack growth step 9 


Table 5 22 Stress ml 

:ensitv factors computed at the end of step 9 

Crack 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Kj m psi^/m 

Kn ps%\/m 

1 

2572 

1 

2 

3 30572489 ' 

-0 402490819 

2 97397514 

-0 428100794 

2 

1 

2572 

1 

2 

2 98920629 

-0 354454145 

3 22757431 

-0 445701305 
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step 10 



Figure 5 21 Problem 6 crack growth step 10 


Table 5 23 Stress intensity factors computed at the end of step 10 


Crack 

No of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Ki in psi\/m 

Kjr psiVm 

1 

2560 

1 

2 

3 65941563 

-0 513562382 

3 21282555 

-0 270626624 

2 

2560 

1 

2 

3 17263255 

-0 394161767 

3 71137782 

-0 449254288 


Dtscusston 

1 A close look at the problem reveals that the crack positions are unaltered by a rotation of 
the domain by 180° So one would expect the cracks to propagate m a self similar manner 
The numerical results demonstrate that indeed it is true 
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2 As the cracks grow, the interaction between the two cracks increases leading to an increase 
in the observed stress intensity factors 

3 The cracks tend to turn towards each other and towards the fixed boundarj- points at the 
two ends This confirms with the observed photoelastic data for similar experiment as well 
as to the intuitive reasoning that near the boundary, the crack will tend to grow towards 
the region of high stress concentration 

4 The mixed mode stress intensity factor is mainly ATj dominated Mode II is important 
from the point of view of the orientation of the crack growth direction The crack turns 
primarily due to Ku 
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Chapter 6 
Closure 


6.1 Summary 

In the present work, we developed a general computational code to model crack growth The 
following conclusions can be drawn based on the present work 

1 The presence of multiple singularities m the domain has been successfully addressed How- 
ever more detailed analysis is to be done to improve the results further 

2 The generalized finite element is superior than the standard finite element solution 

3 The present approach has been extended to arbitrary shaped crack 

4 Mesh refinement is necessary in some cases, as in problem 1, where due to small crack 
length the local disturbance from the nearby cracktip is affecting the solution 

5 The stress intensity factors calculated from the interaction integral computed along the 
two layer of elements around the crack tip yield better results 

6 Higher order shape functions along with little mesh refinement could be very effective tool 
to characterise the crack parameters 


6.2 Suggestions for future work 

1 The major difficulty in crack growth modelling in any of the finite element approach is 
to generate a mesh at each step conforming to the crack surface The present study 
uses automatic mesh generator in crude form The mesh generator should be improved by 
incorporating sub-domain wise meshing The whole domain could be divided into a number 
of sub domains such that the crack surface falls on one of the boundary This excludes the 
possibhty of any internal cut-outs 
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2 In the present work, we have enriched only the node sitting at the crack tip This could 
be extended to enrich the whole zone around the crack tip This will improve the results 
immensely 

3 The enrichment functions used in the present method is as proposed by Fleming et al [4] 
However, much better results could be obtained by using the exact functions as derived for 
re-entrant corner 

4 Mesh refinement might be necessary m some cases particularly when the cracks come 
closer to each other This would improve the results immensely Instead of mesh refinement 
throughout the whole domain, it is better to refine only those elements which are associated 
with the node sitting on the crack tip 

5 In the present study contour form of interaction integral is applied to extract stress intensity 
factor m mixed mode analysis However, domain form of interaction integral based on cut- 
oflf function method (see [12]) or contour form of interaction integral around an annular 
ring surrounding the crack tip (see [12]) would yield better results 
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